Mapping the energy surface of PbTiO,3 in multidimensional electric-displacement 

space 
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In recent years, methods have been developed that allow first-principles electronic-structure cal- 
culations to be carried out under conditions of fixed electric field. For some purposes, however, it 
is more convenient to work at fixed electric displacement field. Initial implementations of the fixed- 
displacement-field approach have been limited to constraining the field along one spatial dimension 
only. Here, we generalize this approach to treat the full three-dimensional displacement field as 
a constraint, and compute the internal-energy landscape as a function of this multidimensional 
displacement-field vector. Using PbTiOa as a prototypical system, we identify stable or metastable 
tetragonal, orthorhombic and rhombohedral structures as the displacement field evolves along [001], 
[110] and [111] directions, respectively. The energy minimum along [001] is found to be deeper than 
that along [110] or [111], as expected for a system having a tetragonal ground state. The barriers 
connecting these minima are found to be quite small, consistent with the current understanding 
that the large piezoelectric effects in PbTiC>3 arise from the easy rotation of the polarization vector. 

PACS numbers: 77.80.-e,71.15.-m 



I. INTRODUCTION 

Since their introduction almost a decade agoji^ meth- 
ods for carrying out first-principles electronic-structure 
calculations under conditions of fixed electric field £ have 
found wide application in the study of the dielectric, 
piezoelectric, and ferroelectric behavior of materials^— 
More recently, variants of this approach, in which the 
electric polarization^ or the electric displacement field 8 " 
is taken as the fundamental variable instead, have been 
introduced. Fixing the displacement field D has the 
intuitive interpretation of imposing open-circuit elec- 
trical boundary conditions, which is often especially 
useful for studying layered geometries such as metal- 
oxide interfaces*^ and ferroelectric capacitors-^ and 
superlatticesj 12 : 13 

There are several reasons why the fixed-D approach 
is advantageous for such applications. In a superlatticc 
structure, the local polarization P and electric field £ 
can vary from one layer to another, so their overall spa- 
tial average is not a fundamental quantity. In contrast, 
D is constant throughout the system, since free charge is 
assumed to be absent. (Here, £, P and D refer to the field 
components in the stacking direction.) Therefore, choos- 
ing D as the fundamental electrical variable is especially 
practical because it makes it possible to decompose the 
equation of state of a layered structure into the sum of 
contributions from the individual building blocks— with 
each of these contributions depending only on the local 
environment. Moreover, long- and short-range electro- 
static interactions are separated efficiently. Such an ap- 
proach facilitates the design of superlattice structures of 
desired functionality by appropriate choice of the stack- 
ing sequence i 12 ' 13 

Even for bulk crystals, it turns out that working at 
fixed D is more convenient than working at fixed £ when 
considering materials that have ferroelectric instabilities. 



(We now refer to three-dimensional field vectors.) The 
reason is that, in the region of the energy landscape near 
the unstable paraelectric configuration, the system is un- 
stable to transformation into one of the ferroelectric do- 
main states when working at fixed £. The unstable para- 
electric region of the E(P) energy landscape is thus inac- 
cessible using this approach. When working at fixed D, 
on the other hand, the internal energy U(D) is typically 
found to be a single- valued function of D,— thus allow- 
ing access to the entire electric equation of stated More- 
over, the second derivative of U with respect to D is di- 
rectly related to the inverse capacitance of the material. 8 
When this quantity goes negative, it is a signature of the 
appearance of a ferroelectric instability, and indeed the 
magnitude of its negative value has been shown to be 
an insightful indicator of the strength of the ferroelec- 
tric instability* 8 -^ For example, it can play an important 
role in determining the critical thickness for ferroelec- 
tricity in capacitor nanostructurespAi and its dependence 
on material structure and composition can be helpful in 
understanding the origin of the ferroelectric instability. 

Up to now, applications of the fixed-displacemcnt- 
field approach have been carried out using the private 
LAUTREC code package, 16 in which the fixed- D con- 
straint can be applied in only one Cartesian direction. 
In the present work, we have implemented the multidi- 
mensional fixed-D method in the context of the open- 
source ABINIT code package^ and demonstrates that it 
can be used to calculate the energy surface throughout 
the region of instability associated with the ferroelectric 
behavior, using PbTiC>3 as our prototypical material. By 
mapping the internal-energy landscape in the full three- 
dimensional D space, we can easily compare the ener- 
gies of the competing ferroelectric states, trace the paths 
connecting these states, and compute the energy barri- 
ers along these paths. This approach therefore gives us 
a powerful tool for characterizing the ferroelectric be- 
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havior of a given material in all of its three-dimensional 
complexity in D (or £ or P) space. 

The paper is organized as follows. In Sec. [ill we briefly 
review the formalism for the fixed-D calculations and 
describe our implementation. In Sec. IIII1 we first test 
the implementation and compare with a previous cal- 
culation on PbTi03. We also explore the internal en- 
ergy landscape of PbTiC>3 in multidimensional D space, 
and map out the relationships between the various field 
variables and their corresponding energy functionals. Fi- 
nally, Sec. IIVI contains a summary and conclusions. 



II. FORMALISM AND METHODOLOGY 

We begin by briefly reviewing the fixed-5 and fixed-D 
formalisms. For the former;^ the natural energy func- 
tional is the electric enthalpy 

F(£,v) = E KS (v)-Sl£--P(v), (1) 

while the fixed-D method 8 is naturally formulated in 
terms of the internal energy 

U(B,v)=E KS (v) + ^{r>-4irP(v)} 2 . (2) 
sir 

In these equations Q is the cell volume, v denotes the 
internal (ionic and electronic) coordinates, and -Eks is 
the ordinary zero-field Kohn-Sham energy functional. 

To implement the fixed-D method based on the above 
formalism, we have modified the open-source ABINIT 
code package, in which the fixed-£ calculation is already 
available*!*^ The key step is to update the electric field 
£ after each SCF iteration according to 

£ n+1 = A(D - 47rP„) + (1 - X)£„, (3) 

where P„ and £ n are polarization and electric field af- 
ter the n'th SCF iteration and A is a damping parame- 
ter used to control the convergence speed. The iteration 
continues until the normal SCF convergence criterion is 
reached, and in addition, |D — 47rP„ — £ n \ becomes less 
than a given tolerance. 

Our calculations were performed within density- 
functional theory in the local-density approximation 19 
using norm-conserving pseudopotentials^ and a plane- 
wave cutoff of 60 Ha. The pseudopotentials were the 
same as those of Ref H. A 6 x 6 x 6 Monkhorst-Pack 
grid^i was used to sample the Brillouin zone. The 
atomic coordinates of the five-atom unit cell were re- 
laxed until all atomic force components were smaller than 
1CP 5 Ha/Bohr, and the cell size and shape was varied un- 
til all stress components were below 10~ 7 Ha/Bohr 3 . 
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FIG. 1: (Color online) (a) Reduced electric field e, and (b) 
internal energy U, as a function of reduced electric displace- 
ment field d along the [0011 direction. The LAUTREC results 
are reproduced from Ref. [g. The solid curve in (a) is a cubic- 
spline fit to the ABINIT results; the solid curve in (b) is the 
numerical integral of the curve in (a). 

III. RESULTS 
A. Testing consistency 

In order to test our implementation and compare the 
results with the previous calculation of Ref. H, we have 
computed the internal energy U and the electric field £ 
as a function of D for the case that D lies along the 
[001] axis. We started the calculation from the relaxed 
cubic structure (lattice constant of 3.88 A) at D — 0, and 
increased D in steps of 0.02 a.u. At each D, the structure 
was fully relaxed with respect to both ionic positions and 
lattice parameters. 

The results are plotted in Fig. [3(a), except that "re- 
duced" field variables are used. That is, for each data 
point we computed the reduced electric field e = c£ and 
the reduced displacement field d — a 2 D/Air, where a and 
c are the x-y and z lattice parameters respectively at 
the given value of D along [001]. In Fig. [U(b), we plot 
the internal energy U vs. reduced displacement field d. 
The previous results of Ref. 8, which were presented in 
terms of the same reduced variables, are also included 
in the plots. It is evident that the agreement between 
our Abinit calculation and the previous Lautrec one is 
excellent. 

An advantage of using reduced variables in Fig. Q] is 
that d, e, and U are related by the exact relation 



even when the relaxation of the lattice vectors with D is 
included. 8 (By contrast, the corresponding relation be- 
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FIG. 2: (Color online) Internal energy surface U(D) of 
PbTiOa plotted for D lying in the plane spanned by the [100] 
and [010] directions, (a) Perspective plot, (b) Contour plot, 
with 3.0 meV separating contour levels. The minimum at 
T and saddle point at O represent spontaneously polarized 
tetragonal and orthorhombic states respectively. 



FIG. 3: (Color online) Internal energy surface U(D) of 
PbTiOs plotted for D lying in the plane spanned by the [110] 
and [001] directions, (a) Perspective plot, (b) Contour plot, 
with 3.0 meV separating contour levels. The minima at T 
and O and the saddle point at R represent spontaneously po- 
larized tetragonal, orthorhombic, and rhombohedral states, 
respectively. 



tween D, £ , and U would require correction terms arising 
from derivatives of lattice vectors and cell volume.) To 
test whether our data is consistent with Eq. Q, we car- 
ried out a cubic spline fit of the s vs. d results to produce 
the solid curve shown in Fig.QJa), and then integrated it 
numerically to obtain the solid curve shown in Fig. [T}3. 
It can be seen that the internal energies coming directly 
from the Abinit calculations coincide quite precisely with 
those predicted by the application of Eq. 

These tests, which were repeated along other directions 
such as [110] with equally good results, confirm the cor- 
rectness and internal consistency of our computational 
implementation. 



B. Internal energy in multidimensional D space 

We now explore the internal energy landscape of 
PbTiC>3 in three-dimensional D space. While our calcu- 
lations can be done for arbitrary D, we restrict ourselves 
for presentation purposes to two-dimensional planes in D 
space. We begin with the case of D lying in the x-y plane. 



Recall that we already obtained relaxed solutions for a se- 
ries of -D[ioo] values increasing in increments of 0.02 a. u., 
up to 0.2 a. u., as described in the previous section. For 
each of these values of -D[ 10 o], we use this solution as a 
starting structure as we apply -D[ 010 ] in steps of 0.02 a.u., 
up to 0.2 a.u., while keeping Z?[ioo] constant. In other 
words, each equilibrium state is used as starting guess 
for its neighbor along y. In this way U(T>) is obtained 
on the specified grid in the quadrant with £)[ioo] > and 
-D[oio] > 0, and then symmetry is used to obtain the full 
internal-energy landscape. 

The result is plotted in Fig. [5] Stationary points in 
such a diagram correspond to states with £ = 0. The four 
minima in the ±[001] and ±[010] directions correspond to 
four of the six tetragonal (T) ground states, at which P 
takes on its spontaneous value P s and D = 4irP s . There 
are also four saddle points along the [110] and related 
directions, corresponding to states of orthorhombic (O) 
symmetry. The fact that these are higher in energy than 
the tetragonal minima just reflects the well-known fact 
that PbTiC>3 has a tetragonal ground state at T = 0. 
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FIG. 4: (Color online) Dark blue: isosurfaces of U(D) for 
PbTiOa plotted at U = — 44meV. Contoured plane: alterna- 
tive view of data of Fig. EJb). 



The figure also shows that the energy differences are 
quite small along the valley connecting states T — > O 
— > T, relative to the large barrier that would have to be 
overcome if one were to pass through the origin of the 
figure. This is consistent with the work of Cohen and 
collaborators^— who pointed out how the easy rotation 
of the orientation of the polarization could lead to large 
piezoelectric responses, even if the energy landscape is 
relatively stiff with respect to changes in the magnitude 
of the polarization. Regarding the piezoelectric response 
of PbTiC>3 starting from its tetragonal ground state, the 
relative flatness of the energy landscape near state T in 
Fig. [2] along that path that would lead from T -> O ->• T 
is heuristically consistent with a relatively large observed 
value of the eis piezoelectric constant in this material.™ 

Fig. [3] shows corresponding plots for the plane span- 
ning the [110] and [001] directions in three-dimensional 
D space. The minima at the top and bottom of panel (b) 
are tetragonal states equivalent to those in Fig. The 
apparent local minima at left and right in panel (b) are 
actually saddle points in three-dimensional D space, and 
correspond to the orthorhombic states already discussed 
in connection with Fig. [3] We now also see four equiva- 
lent saddle points corresponding to a rhombohedral (R) 
state with the polarization along [111] or related direc- 
tions; these are points of double instability, in the sense 
that the Hessian of U(D) has two negative eigenvalues. 
Once again, it is evident that there is a polarization ro- 
tation path via T — > R — > O that is relatively low in 
energy compared to a direct polarization reversal pass- 
ing through the paraelectric maximum at the origin. 

The properties of PbTiC>3 at the tetragonal, or- 
thorhombic and rhombohedral phases are summarized in 
Table |H The tetragonal phase has the lowest internal en- 
ergy, followed by the orthorhombic and then the rhombo- 
hedral phase, consistent with previous studies and with 
the well-known fact that the ground state of PbTiC>3 is 
tetragonal. It can also be seen that the displacement field 
D min minimizing U and the corresponding spontaneous 
polarization P s = D m i n /47r decrease when going from T 



TABLE I: Properties of PbTi0 3 in tetragonal (T), orthorhom- 
bic (O), and rhombohedral (R) phases. The cubic phase is 
chosen to define the zero of the internal energy U. -D m i n is 
the displacement field at which U is a minimum, and P s is the 
corresponding spontaneous polarization. The lattice vectors 
are also given. 
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a, b, c 




(meV/cell) 


(a.u.) 


(C/m 2 ) 


(A) 


T 


-47.78 


0.17 


0.78 


a =3.85, c =4.03 





-39.80 


0.15 


0.68 


a = b =3.92, c =3.86 


R 


-37.23 


0.14 


0.65 


a=3.90, a=89.63° 



to O to R, as might have been guessed from the fact that 
T is the deepest minimum while R is the shallowest. 

To demonstrate the full 3D capabilities of the method, 
we also compute U(D) on a 3D mesh of D values in in- 
crements of AD[ioo] = AD[ io] = AD[ooi] = 0.04 a.u. 
and use this data to plot the internal-energy isosurface 
at U = — 44meV shown in Fig. |4j The six isosurfaces 
surround the six equivalent internal-energy minima (at 
U = -47.3 meV) located along the [001], [010] and [001] 
axes in D space, providing a clear visualization of the six 
possible spontaneously-polarized domain states in tetrag- 
onal PbTiC>3. The contour plane at -D[ooi] = shows 
similar information as in Fig. EJ^b). These results demon- 
strate the robustness of our implementation even in a 
fully multidimensional context. 



C. Relations between field variables 

The electric equation of state of a given crystalline in- 
sulator is given by specifying the relation between any 
two of the three field variables P, D, and £, as for exam- 
ple by the functions £ (D), D(P), or P(£ ). It is straight- 
forward to convert between these using D = £ + 4-7rP. In 
our approach we obtain £(D) directly, and then generate 
D(P) or P(£ ) by numerical manipulation. In general the 
electric equation of state is a vector function of a vector, 
so to simplify our presentation we have calculated and 
plotted the electric equations of state only for cases in 
which all the field variables are constrained to lie along 
either the [001], [110], or [111] axis. 

The results are presented in Fig. [5J with the electric 
equations of state of the form £{D), D(P), and P{£) 
plotted in panels (a-c) respectively. The different curves 
correspond to plots along the [001], [110], or [111] axis. 
We have also marked several special states on the dia- 
grams for the case of the fields being along the [001] axis. 
Proceeding from 1 — > 2 — > 3 (or equivalently by symme- 
try from 1 — > 4 — > 5), we pass from the cubic paraelectric 
State 1 to the spontaneously polarized tetragonal ferro- 
electric State 3. As the displacement field D increases 
in panel (b), the polarization P increases nearly linearly, 
while the electric field Z in panel (a) at first decreases 
to its minimum value (which defines State 2) and then 
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FIG. 5: (Color online) Electric equations of state of the form £(D) (a), D(P) (b), and P(£) (c), plotted for fields constrained 
to lie along the [001], [110] or [111] directions. All units area.u. Numbered dots on the [001] curves indicate special states as 
described in the text. 



increases and passes through zero at State 3. The plot 
of P(£ ) in panel (c) takes the form of a hysteresis curve, 
but the portions of this curve in the region 4 — > 1 — > 
2 are unstable and therefore inaccessible under fixed-£ 
boundary conditions. This is a region in which the di- 
electric permittivity x = dP/d£ is negative. The ability 
of the fixed- D method to explore this region of instabil- 
ity, which cannot be done using the fixed-£ method, is 
one of the important advantages of working at fixed 

Similar behaviors appear for D along the [110] or 
[111] directions. For these directions, States 1-5 are not 
marked, but are defined in the same way. The region 
of instability from State 1 to 2 is almost the same along 
all three directions, as can be seen in panels (a) and (c). 
In all cases, D is very nearly linear in P in panel (b). 
The main difference comes in the D and P values at the 
spontaneously-polarized State 3, which, as already dis- 
cussed in connection with Table HI increase as we go from 
the [111] (R) to the [110] (O) to the [001] (T) directions. 
Also, £ increases much faster along [110] and [111] than 
along [001] after crossing the unstable State 2 in panel 
(a). 

We emphasize that the hysteresis curves shown in 
Fig. [SJc) should not be compared directly with experi- 
ment. They correspond to the theoretical intrinsic hys- 
teretic behavior that would occur if the entire crystal 
would switch coherently (i.e., maintaining full crystal pe- 
riodicity at all times) from +P S to — P s on a path cross- 
ing through P = 0. This is a highly unrealistic picture 
of real ferroelectric switching, which usually proceeds by 
the motion of a domain wall between domains of differ- 
ent orientation of the polarization. 25 Our intrinsic coer- 
cive field of 2.5MV/cm for the [001] case, which can be 
obtained from £ at State 2 in panel (a) or (c), is sure 
to be very much larger than the experimental one, which 
is most likely determined by pinning of domain walls to 
defects or by nucleation phenomena. 



D. Relations between energy functionals 

It is also instructive to see how the energy function- 
als behave as one traverses the trajectories shown in 
Fig. [5] Recall that the energy functionals that are nat- 
urally associated with field variables D, P, and £ are 
U(D), Eks(P), and F{£)\ these are plotted in panels 
(a-c) of Fig. [51 respectively. Each plot again shows the 
behavior for fields constrained along either [001], [110], 
or [111], and the special States 1-5 are again marked for 
the [001] case. 

The U{D) and Ek$(P) plots in panels (a-b) look re- 
markably similar after a rescaling of the horizontal axis 
by a factor of Air. This is not surprising, since States 
1, 3, and 5 have £ = 0, and are thus guaranteed to be 
extrema and to appear in exactly the same place in both 
panels (after 47r rescaling). The inflection points corre- 
sponding to States 2 and 4 are not located in quite the 
same place in both diagrams, but it is difficult to tell this 
by eye. Both panels clearly show a qualitatively similar 
double-well potential. 

The T{£ ) curve in Fig. [HJc) looks complex, but its 
consistency with Fig. [5jc) can be checked through the 
relation P = —dF/d£, which follows from Eq. ([1]). Go- 
ing through the unstable region from State 1 to State 2 
in Fig.lHJc), £ becomes negative while —dJ-/d£ becomes 
positive, consistent with the corresponding behaviors of £ 
and P in Fig. [5jc) . Then in the metastable region from 
State 2 to State 3, and into the stable region beyond 
State 3, £ returns to zero and then goes negative while 
—8J r /d£ continues to grow more positive, again consis- 
tent with Fig. [5][c). The fact that T diverges to — oo as 
\£\ becomes large may look strange, but it is the normal 
behavior of T{£). For example, for a simple linear di- 
electric, U(D) and Eks(P) are simple upright parabolas, 
while J-(£ ) is an inverted parabola^ and the asymptotic 
behaviors at large |£| are similar to those appearing in 
Fig. El 
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FIG. 6: (Color online) Energy functionals vs. their respective natural variables along [001], [110] and [111] directions: U(D) 
(a), Eks(P) (b), and J-(£) (c). All field variables are in units ofa.u. Dots labeled 1-5 indicate the same special states as in 
Fig.! 



IV. SUMMARY AND CONCLUSIONS 

In summary, we have demonstrated the possibility of 
carrying out first-principles density-functional calcula- 
tions under boundary conditions in which all three com- 
ponents of the electric displacement field D are fixed. 
We have implemented the method in the open-source 
ABINIT software package. 

Using PbTi03 as a prototypical system, we have ex- 
plored the internal-energy landscape as a function of the 
full three-dimensional displacement-field vector. We have 
identified the minimum-energy tetragonal, orthorhombic 
and rhombohedral structures, and confirmed that the 
computed properties agree with previous first-principles 
studies. Our results allow for easy visualization of the 
low-energy paths for polarization rotation, known to be 
associated with large piezoelectric responses in this class 
of compounds. We have also presented the electric equa- 



tions of state relating £, D and P, as well as the corre- 
sponding energy functionals, along symmetry lines in D 
space. 

We hope that this fixed-D approach may be useful for 
exploring the internal energy landscape of more compli- 
cated ferroelectric and dielectric materials, in bulk or 
superlattice form, as well as for studying domain-wall 
properties, piezoelectric and flexoelectric responses, field- 
driven phase transitions, and other phenomena in this 
important class of materials. 
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